install.packages("survival"); install.packages("coxme")

library(survival); library(coxme)




#setwd("~/Dropbox/Job Market Paper/IO Revisions/Replication_Files/Replication_Main_Text/Figure_3_&_4")


dur_dat<-read.csv(file="figures_3_4.csv")



###### Figure 4 #########



Surv_1<-Surv(time=dur_dat$start,time2=dur_dat$stop,event=dur_dat$fail_1)

fit.1<-coxph(Surv_1~logArea,data=dur_dat)


new_data_hi<-data.frame(quantile(dur_dat$logArea,.75))

names(new_data_hi)<- "logArea"

new_data_lo<-data.frame(quantile(dur_dat$logArea,.25))

names(new_data_lo)<-"logArea"




pdf(file="Figure_4.pdf")

plot(survfit(fit.1,newdata=new_data_hi,se.fit=F),mark.time=FALSE,col="dodgerblue4",lwd=2,main="",xlab="Time",ylab="Survival Function",axes=F)


times_1<-c(survfit(fit.1,newdata=new_data_hi,se.fit=T)$time)

times_2<-sort(times_1,decreasing=TRUE)

fails_1<-c(survfit(fit.1,newdata=new_data_hi,se.fit=T)$upper)

fails_2<-c(survfit(fit.1,newdata=new_data_hi,se.fit=T)$lower)


lines(stepfun(times_1,c(1,fails_2)),col="dodgerblue4",lwd=1,lty=2,cex.points=0)


lines(stepfun(times_1,c(1,fails_1)),col="dodgerblue4",lwd=1,lty=2,cex.points=0)

#####

times_1a<-c(survfit(fit.1,newdata=new_data_lo,se.fit=T)$time)

times_2a<-sort(times_1,decreasing=TRUE)

fails_1a<-c(survfit(fit.1,newdata=new_data_lo,se.fit=T)$upper)

fails_2a<-c(survfit(fit.1,newdata=new_data_lo,se.fit=T)$lower)


survs<-survfit(fit.1,newdata=new_data_lo,se.fit=F)$surv


lines(stepfun(times_1,c(1,survs)),lwd=2,col="gold3",cex.points=0)


lines(stepfun(times_1a,c(1,fails_2a)),col="gold3",lwd=1,lty=2,cex.points=0)


lines(stepfun(times_1a,c(1,fails_1a)),col="gold3",lwd=1,lty=2,cex.points=0)


axis(2)

axis(1,seq(0,700,by=100))

text(115,.375,".75th",cex=.75)


text(225,.525,".25th",cex=.75)


dev.off()


####################################################################################


fit.2<-coxme(Surv_1~logArea +(1|ID)+(1|Year),data=dur_dat)

fit.3<-coxme(Surv_1~logArea +(1|ID)+(logArea|cent)+(1|Year),data=dur_dat)

fit.4<-coxme(Surv_1~logArea +(logArea|fifty)+(1|ID)+(1|Year),data=dur_dat)


################################
#######Top Left of Fig 3###########
################################


pdf("Figure_3_top_left.pdf")

layout(matrix(c(1,1,3,2),2,2) )


fits<-c(coef(fit.1)[1],coef(fit.2)[1],coef(fit.3)[1],coef(fit.4)[1])

se<-sqrt(c(vcov(fit.1),vcov(fit.2),vcov(fit.3),vcov(fit.4)))

plot(coef(fit.1),xlim=c(1,6),type="n",xlab="",ylim=c(-0.005,.275),axes=F,ylab="Estimate Size")

grid()

text(1,.25,expression(hat(mu[gamma])), cex=1.5)

for(i in 1:4){

points(i+1,fits[i],pch=19,cex=1.25)

lines(c(i+1,i+1),c(fits[i]+1.96*se[i],fits[i]-1.96*se[i]),lwd=2)

lines(c(i+1-.05,i+1+.05),c(fits[i]+1.96*se[i],fits[i]+1.96*se[i]),lwd=2)

lines(c(i+1-.05,i+1+.05),c(fits[i]-1.96*se[i],fits[i]-1.96*se[i]),lwd=2)

}

axis(2,seq(-.05,.3,by=.05))

abline(h=0,lty=2)

axis(1,seq(1,6),labels=c("","Model 1","Model 2","Model 3","Model 4",""))
axis(1,seq(1,6),labels=c("","","","Cox Mixed Effects","",""),tick=F,line=.95)
axis(1,seq(1,6),labels=c("","","","","Time-Varying Slopes",""),tick=F,line=1.75)



dev.off()












###################################
#######Bottom Right of Fig 3###########
###################################

pdf("Figure_3_bottom_right.pdf",width=7, height=5)

Year2<-seq(1150,1800,50)

plot(Year2,ranef(fit.4)$fifty,ylim=c(-.1,.12),pch=19,axes=F,xlab="Time-Varying Parameter",ylab="Estimate Size",main="Model 4")

grid()


ci_4_up<-sqrt(fit.4$vcoef$fifty)*1.96+ranef(fit.4)$fifty
ci_4_lo<- -1*sqrt(fit.4$vcoef$fifty)*1.96+ranef(fit.4)$fifty

for(i in 1:14){
	
lines(c(Year2[i],Year2[i]), c(ci_4_lo[i],ci_4_up[i]))

lines(c(Year2[i]-5,Year2[i]+5), c(ci_4_lo[i],ci_4_lo[i]))

lines(c(Year2[i]-5,Year2[i]+5), c(ci_4_up[i],ci_4_up[i]))

}


abline(h=0,lty=2)


axis(2,at=seq(-.25,.25,by=.05))

axis(1,at=Year2,labels=c(expression(hat(gamma)[1150]),expression(hat(gamma)[1200]),expression(hat(gamma)[1250]),expression(hat(gamma)[1300]),expression(hat(gamma)[1350]),expression(hat(gamma)[1400]),expression(hat(gamma)[1450]),expression(hat(gamma)[1500]),expression(hat(gamma)[1550]),expression(hat(gamma)[1600]),expression(hat(gamma)[1650]),expression(hat(gamma)[1700]),expression(hat(gamma)[1750]),expression(hat(gamma)[1800])))



dev.off()




################################
#######Top Right of Fig 3###########
################################



pdf("Figure_3_top_right.pdf",width=7, height=5)




Year1<-seq(1200,1800,by=100)

plot(Year1,ranef(fit.3)$cent,ylim=c(-.10,.11),pch=19,axes=F,xlab="",ylab="Estimate Size",main="Model 3")

grid()


ci_3_up<-sqrt(fit.3$vcoef$cent)*1.96+ranef(fit.3)$cent
ci_3_lo<- -1*sqrt(fit.3$vcoef$cent)*1.96+ranef(fit.3)$cent

for(i in 1:7){

lines(c(Year1[i],Year1[i]),c(ci_3_lo[i],ci_3_up[i]),lty=1,lwd=2)

lines(c(Year1[i]-5,Year1[i]+5), c(ci_3_lo[i],ci_3_lo[i]))

lines(c(Year1[i]-5,Year1[i]+5), c(ci_3_up[i],ci_3_up[i]))

}

abline(h=0,lty=2)


axis(1,at=Year1,labels=c(expression(hat(gamma)[1200]),expression(hat(gamma)[1300]),expression(hat(gamma)[1400]),expression(hat(gamma)[1500]),expression(hat(gamma)[1600]),expression(hat(gamma)[1700]),expression(hat(gamma)[1800])))

axis(2,at=seq(-.25,.25,by=.05))


dev.off()



save.image(file="figs_3_4.Rdata")




